/* 
 * File:   CompressedRowStorage.h
 * Author: Hugo
 *
 * Created on October 28, 2013, 12:34 PM
 */

#ifndef COMPRESSEDROWSTORAGE_H
#define COMPRESSEDROWSTORAGE_H

#include "vector.h"
#include "matrix.h"
#include <map>

template <class T>
class CompressedColumnStorage: public Matrix<T>
{
template <class U>    
friend Matrix<U> operator* (const Matrix<U>&,const CompressedColumnStorage<U>&);
template <class U>    
friend const Vector<U> operator*(const Vector<U>&, const CompressedColumnStorage<U>&);

public: 
    
    CompressedColumnStorage(const Matrix<T>& _M);
    CompressedColumnStorage(unsigned int m_row, unsigned int m_col = 0);
    CompressedColumnStorage(const CompressedColumnStorage<T>& _CS);
    
    //Destructor:
    ~CompressedColumnStorage();
    
    //Utility Methods:
    void zeroRowColumn(unsigned int _k);
    double getElement(int i, int j);
    void print() const;
    Matrix<T> operator*(const Matrix<T>&) const;	
    Vector<T> operator*(const Vector<T>&) const;
    Matrix<T> operator + (const Matrix<T>& m) const;
    Matrix<T> operator - (const Matrix<T>& m) const;
    Matrix<T> operator * (const T& a) const;
    Matrix<T> toMatrix();
    Vector<T> isolve( const Vector<T>& _V );
    //TODO: product between two compressed-columns!
    
    //Equality Operators:
    bool operator == (const Matrix<T>& ) const;
    bool operator != (const Matrix<T>& ) const;
    T& operator () (unsigned int row, unsigned int col);
    Vector<T> PCGsolve( const Vector<T>& _V,double error = 0.0001);
    Vector<T> CGsolve( const Vector<T>& _V, double e);
    void initData();
    //TODO: get subscript operator working!
private:

	int m_dim;
        bool isFactorized;
	IVector col_ptr;
	IVector row_index;
	Vector<T> value;
	Matrix<T> pre_multiply(const Matrix<T>&)const;
	Matrix<T> post_multiply(const Matrix<T>&) const;
	Vector<T> pre_multiply(const Vector<T>&)const;
	Vector<T> post_multiply(const Vector<T>&) const;
        void ichol();
        map<int,T> *col_map;
        
        
};




template <class T> inline
CompressedColumnStorage<T>::CompressedColumnStorage(const Matrix<T>& _M)
: Matrix<T>(_M),
isFactorized(false)
{
    
    m_dim = _M.getCols();
    int non_zeros = m_dim*(m_dim+1)/2;
    
    
    //total of non_zeros elements:
    for(int j=0; j<m_dim; j++){
        for(int i=0; i<j; i++){
            if(_M(i,j) == 0){
                non_zeros--;
            }
        }
    }
    
    
    row_index.resize(non_zeros);
    col_ptr.resize(m_dim + 1);
    value.resize(non_zeros);
    
    //storage method:
    int index=0;
    col_ptr(0) = 0;
    for(int j=0; j<m_dim; j++){
        for(int i=0; i<=j; i++){
            if(_M(i,j) != 0){
                value(index) = _M(i,j);
                row_index(index++) = i;
            }
            
        }
        col_ptr(j+1) = index; 
    }
    
}

//-----------------------------------------------------

//-----------------------------------------------------

template <class T> inline
CompressedColumnStorage<T>::CompressedColumnStorage(unsigned int m_row, unsigned int m_col):
col_map(new map<int,T>[m_row]),
isFactorized(false),
m_dim(m_row),
col_ptr(m_dim+1)
{

 
}

//-----------------------------------------------------

//-----------------------------------------------------

template <class T> inline
CompressedColumnStorage<T>::CompressedColumnStorage(const CompressedColumnStorage<T>& _CS):
col_map(0),
isFactorized(false),
m_dim(_CS.m_dim),
col_ptr(_CS.col_ptr),
row_index(_CS.row_index),
value(_CS.value)
{

 
}

//-----------------------------------------------------

//-----------------------------------------------------

template <class T> inline
CompressedColumnStorage<T>::~CompressedColumnStorage()
{
 
 
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
double CompressedColumnStorage<T>::getElement(int i, int j){

    // Find where the row starts
    int columnStart = col_ptr(j);

    // Find where the next row starts
    int columnEnd = col_ptr(j+1);

    // Scan the column indices of entries in row i
    for(int k = columnStart; k < columnEnd; k++)
    {
            // if j is the column index of a non-zero, then return it
            if(row_index(k) == i)
                    return value(k);

            // if we have passed j, then entry (i, j) must be a zero
            if(row_index(k) > i)
                    return 0;
    }

    // if we have reached the end of the non-zeroes without
    // find j, it must be the index of a trailing 0
    return 0;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void CompressedColumnStorage<T>::print() const{
    cout << "Values vector:\n";
    for(int i=0; i<value.getLength(); i++){
             cout << value(i) << "\n";
    }
    cout << "\nRow Index Vector:\n";
    for(int i=0; i<row_index.getLength(); i++){
            cout <<row_index(i) << "\n";
    }
    cout << "\nColumn Pointer Vector:\n";
    for(int i=0; i<col_ptr.getLength(); i++){
            cout << col_ptr(i) << "\n";
    }
}
//-----------------------------------------------------

template <class T> inline
void CompressedColumnStorage<T>::ichol() {
    
    Vector<T> iw(m_dim);
    for(int j=1; j < m_dim;j++){
        iw.zero();        
        //Reduce the off-diagonal elements:
           for(int i=col_ptr(j); i<col_ptr(j+1)-1 ; i++){
               iw(row_index(i)) = value(i);
                T sum = 0;
                for(int m = col_ptr(row_index(i)); m< col_ptr(row_index(i)+1) - 1; m++){
                    sum += value(m)*iw(row_index(m));
                }
                value(i) = value(i) - sum;
                iw(row_index(i)) = value(i);
           }
           
           
           //Reduce the diagonal elements
           
           T sum = 0;
           for(int m = col_ptr(j); m < col_ptr(j+1) - 1; m++){
               sum += value(m)*value(m)/value(col_ptr(row_index(m)+1)-1);
           }
           
           value(col_ptr(j+1)-1) = value(col_ptr(j+1)-1) - sum;
           
           
           //Complete reduce:
           
           for(int m = col_ptr(j); m < col_ptr(j+1) - 1 ; m++){
               value(m) = value(m)/value(col_ptr(row_index(m)+1)-1);
           }
           
           
        
    }
    isFactorized = true;
}

//-----------------------------------------------------

template <class T> inline
Vector<T> CompressedColumnStorage<T>::isolve( const Vector<T>& _V )
{
 if(!isFactorized)
     this->ichol();
 Vector<T> r(_V);
 double sum = 0.;
  for(int i=1; i<m_dim; i++){
   sum = 0.;
   for(int j=col_ptr(i); j<col_ptr(i+1)-1; j++){
    sum = sum + value(j)*r(row_index(j));
   }
   r(i) = r(i) - sum;
  }

  for(int i=0; i<m_dim; i++){
   r(i) = r(i)/value(col_ptr(i+1)-1);
  }
 
 
  for(int j=m_dim-1; j>0; j--){
   for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
    r(row_index(i)) -= value(i)*r(j);
   }
  }
  
  return r;
}

//-----------------------------------------------------
template <class T> inline
Vector<T> CompressedColumnStorage<T>::PCGsolve( const Vector<T>& _V, double e)
{
    this->initData();
    CompressedColumnStorage A(*this);
    Vector<T> x(m_dim);
    Vector<T> r(m_dim);
    Vector<T> d(m_dim);
    Vector<T> aux(m_dim);
    Vector<T> dot_aux(m_dim);
    double dot_new;
    double dot_old;
    double alpha;
    double beta;
    int i;
    i = 1;
    x.zero();
    r = _V;
    d = this->isolve(r);
    dot_new = r*d;
    double error = fabs(e*e*dot_new);
    
    
    while(fabs(dot_new) > error & i<10000){
        dot_aux = A*d;
        alpha = dot_new/(d*dot_aux);
        x = x + d*alpha;
        if(i%50 == 0){
        r = _V - (A*x);
        }
        else{
        r = r - (alpha*dot_aux);
        }
        aux = this->isolve(r);
        dot_old = dot_new;
        dot_new = r*aux;
        beta = dot_new/dot_old;
        d = aux + (d*beta);
        i++;
        }
    cout << "Número de iterações:\n" << i;
    return x;
    }
//-----------------------------------------------------
//-----------------------------------------------------
template <class T> inline
Vector<T> CompressedColumnStorage<T>::CGsolve( const Vector<T>& _V, double e)
{
    this->initData();
    Vector<T> x(m_dim);
    x.zero();
    Vector<T> r;
    r = _V;
    Vector<T> d;
    d = r;
    Vector<T> q;
    double dot_new = r*d;
    double dot_old;
    double alpha;
    double beta;
    int i = 1;
    double error = e*e*dot_new;
    
    
    while(dot_new > error & i<10000){
        q = (*this)*d;
        alpha = dot_new/(d*q);
        x = x + d*alpha;
        if(i%50 == 0){
        r = _V - ((*this)*x);
        }
        else{
        r = r - (q*alpha);
        }
        dot_old = dot_new;
        dot_new = r*r;
        beta = dot_new/dot_old;
        d = r + (d*beta);
        i++;
        }
    cout << "Número de iterações:\n" << i;
    return x;
    }
//-----------------------------------------------------
template <class T> inline
Matrix<T> CompressedColumnStorage<T>::pre_multiply(const Matrix<T> &_M) const{
    Matrix<T> result_matrix(m_dim,_M.getCols());
    result_matrix.zero();
    for(int j=0; j<m_dim;j++ ){
        for(int k=0; k<_M.getCols(); k++){
            for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
                result_matrix(j,k) += value(i)*_M(k,row_index(i));
                result_matrix(row_index(i),k) += value(i)*_M(k,j);
            }
            result_matrix(j,k) += value(col_ptr(j+1)-1)*_M(k,j);
        }
    }
    return result_matrix;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T> CompressedColumnStorage<T>::post_multiply(const Matrix<T>& _M) const {
	Matrix<T> result_matrix(_M.getRows(),m_dim);
	result_matrix.zero();
	for(int j=0; j<m_dim;j++ ){
		for(int k=0; k<_M.getRows(); k++){
			for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
				result_matrix(k,row_index(i)) += value(i)*_M(k,j);
				result_matrix(k,j) += value(i)*_M(k,row_index(i));
			}
			result_matrix(k,j) += value(col_ptr(j+1)-1)*_M(k,j);
		}
	}
	return result_matrix;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Vector<T> CompressedColumnStorage<T>::pre_multiply(const Vector<T>& _V) const{
	Vector<T> result_vector(m_dim);
	result_vector.zero();
	for(int j=0; j<m_dim;j++ ){
			for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
				result_vector(j) += value(i)*_V(row_index(i));
				result_vector(row_index(i)) += value(i)*_V(j);
			}
			result_vector(j) += value(col_ptr(j+1)-1)*_V(j);
	}
	return result_vector;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Vector<T> CompressedColumnStorage<T>::post_multiply(const Vector<T>& _V) const{
	Vector<T> result_vector(m_dim);
	result_vector.zero();
	for(int j=0; j<m_dim;j++ ){
		
			for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
				result_vector(row_index(i)) += value(i)*_V(j);
				result_vector(j) += value(i)*_V(row_index(i));
			}
			result_vector(j) += value(col_ptr(j+1)-1)*_V(j);
	}
	return result_vector;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T> CompressedColumnStorage<T>::operator*(const Matrix<T>& _M) const{
	return this->pre_multiply(_M);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Vector<T> CompressedColumnStorage<T>::operator*(const Vector<T>& _V) const{
	return this->pre_multiply(_V);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T> CompressedColumnStorage<T>::operator+(const Matrix<T> &_M) const{
    Matrix<T> result_matrix(m_dim,m_dim);
    result_matrix.zero();
    for(int j=0; j<m_dim;j++ ){
            for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
                result_matrix(row_index(i),j) = value(i) + _M(row_index(i),j);
                result_matrix(j,row_index(i)) = value(i) + _M(j,row_index(i));
            }
            result_matrix(j,j) = value(col_ptr(j+1)-1) + _M(j,j);
    }
    return result_matrix;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T> CompressedColumnStorage<T>::operator-(const Matrix<T> &_M) const{
    Matrix<T> result_matrix(m_dim,m_dim);
    result_matrix.zero();
    for(int j=0; j<m_dim;j++ ){
            for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
                result_matrix(row_index(i),j) = value(i) - _M(row_index(i),j);
                result_matrix(j,row_index(i)) = value(i) - _M(j,row_index(i));
            }
            result_matrix(j,j) = value(col_ptr(j+1)-1) - _M(j,j);
    }
    return result_matrix;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T> CompressedColumnStorage<T>::operator*(const T &_c) const{
    Matrix<T> result_matrix = *this;
    for(int j=0; j<m_dim;j++ ){
            for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
                result_matrix(row_index(i),j) *=_c;
                result_matrix(j,row_index(i)) *=_c;
            }
            result_matrix(j,j) *=_c;
    }
    return result_matrix;
}
//-----------------------------------------------------
template <class T> inline
bool CompressedColumnStorage<T>::operator == (const Matrix<T>& _M) const{
    for(int j=0; j<m_dim;j++ ){
            for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
                if(_M(row_index(i),j) != value(i) || _M(j,row_index(i)) != value(i)){
                    return false;
            }
            }
            if(value(col_ptr(j+1)-1) != _M(j,j)){
                return false;
            }
    }
    return true;
}
//-----------------------------------------------------
template <class T> inline
T& CompressedColumnStorage<T>::operator () (unsigned int row, unsigned int col){
    if(col<row)
        return col_map[row][col];
    else
    return col_map[col][row];
}

//-----------------------------------------------------
template <class T> inline
bool CompressedColumnStorage<T>::operator != (const Matrix<T>& _M) const{
    return !this->operator ==(_M);
}
//-----------------------------------------------------



//-----------------------------------------------------
template <class T> inline
Matrix<T> CompressedColumnStorage<T>::toMatrix(){
	Matrix<T> matrix(m_dim,m_dim);
	matrix.zero();
	for(int j=0; j<m_dim;j++ ){
		for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
			matrix(row_index(i),j) = value(i);
			matrix(j,row_index(i)) = value(i);
		}
		matrix(j,j) = value(col_ptr(j+1)-1);
	}
	return matrix;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void CompressedColumnStorage<T>::initData(){
    unsigned int non_zeros = 0;
    col_ptr(0) = 0;
    for(unsigned int i=0; i<m_dim; i++){
        non_zeros += col_map[i].size();
        col_ptr(i+1) = non_zeros;
    }
    value.resize(non_zeros);
    row_index.resize(non_zeros);
    unsigned int count = 0;
    for(unsigned int i=0; i<m_dim; i++){
    for ( typename std::map<int,T>::const_iterator iter = col_map[i].begin(); iter != col_map[i].end(); ++iter ){
        row_index(count) = iter->first;
        value(count++) = iter->second;
    }
    }
    delete[] (col_map);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void CompressedColumnStorage<T>::zeroRowColumn(unsigned int _k){
    col_map[_k].clear();
    for(unsigned int i=_k+1; i<m_dim; i++){
    col_map[i].erase(_k);
    }
    col_map[_k][_k] = 1;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class U>
const Matrix<U> operator*(const Matrix<U>& _M, const CompressedColumnStorage<U>& _CS) {
	return _CS.post_multiply(_M);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class U>
const Vector<U> operator*(const Vector<U>& _V, const CompressedColumnStorage<U>& _CS) {
	return _CS.post_multiply(_V);
}
//-----------------------------------------------------




typedef CompressedColumnStorage<double>  DCompressedColumnStorage;
typedef CompressedColumnStorage<int>  ICompressedColumnStorage;
        
#endif	/* COMPRESSEDROWSTORAGE_H */
        




